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Abstract. We present the results of our numerical exploration of tlie fractal structure found by 
S.P. Novikov in the problem of the behaviour of magnetoresistance in a normal metal under a 
strong magnetic field. The case we discuss in this paper is the simplest non-trivial one, namely 
the case of 2 Fermi Surfaces that cut the brillouin zone along of the coordinate axes (i.e. Fermi 
surfaces have genus 3). 

1 Introduction 

Remarkable topological properties of the problem of asymptotics of orbits of quasi-momenta in the dual 
lattice of a normal metal under a strong magnetic field have been noticed by S.P. Novikov in early eighties 
[Nov82]. 

After analysis of the system behavior for magnetic fields close to rational [Zor84] and in "generic position" 
[Dyn93a,Dyn97], the following picture has been extracted by S.P. Novikov (see [NM98] for a review and more 
, bibliography) : once a Fermi function (or a set of dispersion relations) has been fixed, on the space of directions 
of the magnetic field is defined a fractal consisting of smooth polygons that generically have a finite number 
of points in common. Every such polygon is labeled by an integer plane (i.e. a Miller index) and to every 
point of it are associated 2 energies. 

The meaning of these data is the following: suppose a metal has Fermi Function e and Fermi Energy 
Ep and we want to know the asymptotic behavior of trajectories of quasi-momenta for some magnetic field 
H, and be / = fc), ei and 62 the Miller index and energies associated to H. Then the answer is that if 
ei < Ep < 62 there are open orbits and they are finite deformation of the straight line of direction H x I, 
while if Ep < ei or 62 < Ep all orbits are closed. 

From these facts it is clear that knowing the zones, their labels and the functions ei and 62 gives us a 
complete knowledge of the asymptotic behavior of trajectories. Even in the most elementary cases though 
it is impossible to get analytical expression for functions and I, so a numerical analysis of the problem is 
necessary. 

When Ep — ei — 62 we can instead get much more complicated "ergodic-like" behavior of trajectories. 
For those directions no label is defined in general, unless they belong to the boundary of some stability zone. 
It is known [DL99] that generically, in the fixed energy picture, the measure of such directions is zero and 
that the set has a "Cantor-like" fractal structure, but nothing is known about their measure in the global 
picture. 
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The following conjecture has been formulated by Novikov: for a generic set of dispersion relations, the 
measure of "ergodic" directions in the global picture is zero and their fractal dimension is between 1 and 2. 

Our numerical study aims to find the first numerical evidence of the existence of this fractal structure 
and to evaluate its Minkowski fractal dimension [Fal97] in the simplest smooth non trivial case, namely the 
function f{x) = cos(a;) + cos(y) + cos(2;). 

We also repeat calculations for a piecewise polynomial function with the same symmetries, which allows 
us to work more from the analytical point of view. 



2 The idea behind the algorithm 

Let us start explaining what we want our algorithm to do. 

Our ingredients are a function /, smooth or piecewise smooth, a non critical value c of the function that 
give rise to the smooth surface = /^^(c) G and a direction (magnetic field) H € MP^ (we disregard 
here any effect concerning magnetic breakdown; we just assume that our magnetic field is strong enough to 
give rise to the phenomenon, i.e. it is at least of the order of 10^ Gauss, and not strong enough to deform 
the Fermi Surface, so that the only free parameter left is its direction). 

The goal is to get the Miller index associated with H, i.e. the homology class of the 2-tori (if they 
exist) on which lie the open orbits corresponding to H (see [NM98] for details). In other words we must find 
somehow three integer numbers I = {l,m,n) G i?2(T^,Z) that represent the integer homology class of some 
2-torus embedded in T^. 

To be able to understand the way to get these numbers, we make the following consideration: assuming 
we fix some rational magnetic field H such that 62 (ff) ^ ei{H), our surface Mg will be split into an even 
number of 2-tori connected through cylinders of closed orbits, a configuration which is homcomorphic to the 
one shown below. In general a surface of genus g will give rise to at most g — I cylinders that in turn will 
separate at most g — 1 ov g — 2 2-tori when g is respectively odd or even. 

It is clear that if we cut the cylinders with any pair of 2-tori with the same homology class of the ones 
on which the open orbits lie we get a bordism between the sets of cycles cut by them on cylinders. In other 
words, such a section would give us a set of 1-cycles, non trivial in M^, so that the sum z of their homology 
classes does not depend on the height at which we choose to intersect. 

The key point for our algorithm is that there is a 1-1 correspondence between this homology class 
z e Hi{M^, Z) and the homology class h G H2{T^, Z) of the 2-tori. 

Let us call i : ^ the embedding of the surface in the 3-torus and i., : Hi{M^,Z) Hi(T^,Z) 
the induced homomorphism on the first homology groups. It is straightforward to observe that z G ker(i*) 
and that the cycles that lie on the 2-tori, i.e. the ones that are sent by i* in the 2-dimensional sublattice of 
Hi{T^,Z) corresponding to h, have intersection number with z and that, on the contrary, all cycles that 
have an intersection number different from with z cannot lie on a single 2-torus. 

It follows hence that the value of h can be found getting the images by of all cycles that have number 
of intersection with z, or equivalently that are symplectically orthogonal to z with respect to the natural 
structure of symplectic vector space of i?i(M^,M). 

For example, let us sec; what happens in the cases we investigated numerically so far: both functions 
we studied, in the range of values (—1, 1), give rise to a genus-3 surface embedded in T"^ with rank 3 and 
their three handles cut symmetrically (see figures 4.1 and 5.1) the six sides of the cube [0, l]'^ (that we use 
as a model for T^). In particular we are going to have just two cylinders that separate two 2-tori, so every 
section by a 2-torus parallel (i.e. with the same homology class) to them will cut exactly one cylinder and 
give rise to exactly one homology class z G ker . 

If we choose as basic cycles the ones coming naturally from the embedding, as shown in picture 5.1, and 

call them and their canonical duals respectively Qx, ciy and and bx, by and bz, we see that sends the 
tti's respectively to (1, 0, 0) , (0, 1, 0) and (0, 0, 1) and sends the 6j's to (0, 0, 0). The cycle z is then a linear 
combination of bi's with integer coefficients, z = lbx + mby + nbz, and the cycles sent in h by i* are the ones 
such that < c,z >= 0, c = I' ax + m'ay + n'az- 



3 The NTC library 
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Figure 2.1 : an example of the splitting in 2-tori and cylinders of a surface of genus 7 under the 
action of a constant magnetic field. 

It is then clear now that the homology class h is represented in H2 {T^ , by the same triple of integers 
that represents z in Hi{AI^,Z), so in this particular case it is enough to find out the intersection number of 
z with the a^'s to get h. 

3 The NTC library 

After the previous discussion it is clear which capabilities we expect from the software we are going to 
use to perform the numerical analysis. It must be able to deal with the topology of a surface, i.e. it must 
have the possibility of dealing with simplexes of dimension 0, 1, 2 and 3, and it must be able to perform 
topological operations like getting a simplicial decomposition of the level set of a function of three variables 
(to get the Fermi surface M^), intersecting two simplicial complexes (to get the 1-dimensional leaves on 
M^), identifying closed curves in (i.e. it must be able to deal with the periodic boundary conditions 
that identify [0, l]'^ with T'^), evaluating intersection numbers between 2 cycles on a surface and finding the 
homology class of loops in T'^. 

When we started working at this project, after thorough search on the InterNet we found several C++ 
libraries able to deal with the topology of 3- and lower-dimensional objects through simplicial decompositions 
("meshes" in computer jargon). None of them of course directly implements the specific functions we needed, 
so we decided to write a C-I-+ library on top of one of the preexisting ones to implement the complex 
topological functions we needed and tried to make the code as much reusable as possible as it seems that 
such a library could be useful in the future also for different numerical topological problems. We called our 
library Novikov Torus Conjecture library (NTC). 

Afte r an accurate examination of all libraries available we chose to use the library VTK (Visualization 



ToolKit, lit tp://www. kitwarc.com ). The main reasons for our choice, aside from the fact that it is free, are 
the availability of its source code, the fast rate at which is developed and the existence of a very active 
mailing list about VTK-related problems and solutions. Moreover, as VTK was intended primarily as a 



visualization tool based on the standard C library OpenGL by SGI ( ^ittp: / /www.opcngl.org ), it easily allows 



us to visualize our surfaces and cycles, making much easier the debugging process. 
Let us now explain in detail how the algorithm works. 

First of all it is important to point out that, as every computer can basically deal just with integer 
numbers, we are able just to explore what happens for 1-rational magnetic fields. This is not a big restriction 
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as the homology class we look for is locally constant with respect to the direction of the magnetic field and 
so, close to every H , there is a rational one that induces the same homology class. 

Let us then fix some 1-rational H . Its critical points generically will be of "figure eight" type, i.e. its 
tails close up to a pair of loops. We call these critical points "fully open" , "half-closed" and "fully-closed" 
using a terminology coming from R'^, or equivalently depending on whether in T'^ both loops are non trivial, 
one is homotopic to or they are both trivial. 

The first thing to do is to check whether it gives rise to the generic situation or not, i.e. whether open 
orbits lie on 2-tori or not. For simplicity we will just analyze the case of genus 3 surfaces, in which just two 
2-tori and two cylinders can appear, that are the only cases studied numerically by now. 

As we said, the 2-tori components are separated by cylinders of closed orbits that have as basis a pair 
of "half-closed" critical saddle points (see for example figure 3.1), i.e. saddles in which two of the four tails 
are connected so that the loop they form is homotopic to in T^. Our strategy then is to find numerically 
the set of all critical points, single out the saddles that play role in the topology of the surface (we disregard 
all saddles that are associated to a center; they are easy to spot as they are half or fully closed and form a 
loop homotopic to in the surface) and count how many of them are half closed. 




c T3). 



If at least three of them, and so all four, are half closed, then we evaluate the homology class in the 
surface of any of the loops and this would give us exactly the homology class of the 2-tori of open orbits. 

If instead one, and so at least two, of them is fully closed or fully open, then there cannot be any rank-2 
2-torus; there could be still open orbits but they will fill rank-1 tori, i.e. cylinders, and these open orbits 
will disappear for any generic small enough perturbation, so we do not register any homology class for these 
directions. 

This algorithm is implemented through the construction of several classes. The main classes are ntcFoli- 
ation and its subclass ntcPlaneFoliation, that contain all parameters (Fermi function, energy, magnetic field 
and resolution of samplings) and functions able to produce and classify the critical points of the foliation and 
to get the critical leaf (our terminology comes from the universal covering R'^ : in T"^ for 1-rational magnetic 
fields critical leaves will be always fully closed, but the ones leading to open orbits will have just one of the 
loops homotopic to in T"^ so in the covering such saddle is half-open). 

For several reasons we implemented two different ways to get the critical leaves. One works by obtaining 
the level lines of the Fermi function restricted to the plane perpendicular to H passing through the singular 
point. To get the full picture of the intersection we restrict our sampling to an opportunely chosen paral- 
lelogram spanned by a 1? basis of the 2-dimensional lattice given by the intersection of 1? with the plane 
perpendicular to H , so that we get a picture that glues nicely on the boundary. 



4 Study of the trigonometric function 
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This procedure tends to need too much RAM when the components of H get big, say around 400, 
because the area of the basic paraUelogram tends to get too big. We implemented therefore a second way 
that needs just the simphcial decomposition of the Fermi surface in the cube [0, 1]'^. The critical leaf now is 
obtained starting cutting by the plane passing through the critical point and then following the loop. When 
the trajectory reaches the boundary of the cube the coordinates of the equivalent point are evaluated and a 
new plane is taken passing through that point. The process stops when the loop comes back for the second 
time to the critical point (as there are two critical loops for every critical point). 




Figure 3.2 : an example of the level set of the Fermi function obtained through the second kind of 
algorithm implemented. Just the critical loop homotopic to (in T'^) is shown. 

In both cases the critical leaf is at the end copied in an object of the class ntcPrimitiveCell that contains 
methods to deal with loops in T'^. If just one of the two loops is homotopic to zero in T"^, then the methods 
for evaluating the intersection numbers of that cycle with the 3 cycles in the kernel of are called (they 
are contained in the class ntcImplicitFunction that contains all data about the function) and the result gives 
automatically the searched homology class. 

A few other complementary class are also implemented to deal with leaves and functions. The complete 
documentation for the NTC library t ogether with the source code is available at the InterNet address 
http://www.math.umd.edu/~rdl/ntc/ . 

4 Study of the trigonometric function 

The function f{x, y, z) — cos(27ra;) + cos(27r?/) + cos(27rz) is the simplest trigonometric function that 
gives rise to a non trivial (i.e. rank 3) embedding of a surface in T'^ and the only one that had been studied 
so far. 

The only critical values of / are ±3 and ±1, so all level sets — f^^{E) are homeomorphic to spheres 
for E e (—3, —1) U (1, 3). The level set Mq shown in figure 4.1 shows that for E E (—1, 1) all level sets are 
genus-3 surfaces embedded with rank 3 in T"^. 

In particular this means that every generic foliation of Me will have at least four saddle points, and that 
all saddles but four will be associated to some center and hence will be homotopic to in the surface. These 
critical points, that we call "topological" as their origin is due to the topology of the surface and not to the 
particular embedding, are at the base of two cylinders that divide two 2-tori for every generic direction of 
the magnetic field. 

Each level surface of this function is invariant under the symmetry group of the cube; this action in 
turn induces an action on §^ under which the fractal picture is invariant, so it is enough for us to analyze 
its structure in one of the 48 domains in which the action subdivides S^. 

In the projective chart of corresponding to the plane z — 1, one of these domains is the triangle 
2^ < y C [0, 1]^, so we will refer just to the square [0, 1]^ as our "phase space" from now on. On this 
square the picture of stability zones will be symmetric with respect to the diagonal, fact that will be used 
as consistency check of our algorithm. 
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Figure 4.1 : the surface Mq cut by a plane passing through a critical point. 

Another symmetry, due to the fact that the cosine is an even function, impHes that aU level sets 
are symmetric with respect to the origin. It follows that the four topological saddles are divided in two 
symmetrical pairs that define one cylinder each. In fact we can assign to every closed (in M'^) orbit a plus 
or minus sign, according to the fact that it bounds a region where / assumes respectively values smaller or 
bigger than the one assumed on the loop (these two different kinds of loop are called "electrons" and "holes" 
in physics literature). 

This sign is invariant by homotopy, so the same sign is associated to the whole cylinder and is shared 
by the two critical loops at the two bases. As the symmetry respect to the center do not switch this sign, it 
is clear that every pair of symmetric critical points defines one of the two cylinders. 

Finally, the identity cos(27rx) = — cos [27r(l/2 — x)] induces a symmetry between different level surfaces, 
namely the surface is obtained from M_c through a translation and a reflection respect to the origin. As 
the foliation PiH^ — const is invariant by these two operations, it is clear that the existence of open orbits 
at energy c implies the existence of open orbits at energy — c, so that the interval for which any direction 
gives rise to open orbits (that is closed, connected and non empty by [Dyn97]) has the form [—E,E]. 

The surface Mq hence plays a very special role, as at energy c = every direction gives rise to open 
orbits and so every "stability zone" reaches here its biggest size. This means that to study the fractal on 
§^ corresponding to this function is enough to study the level c — 0, while in general it would be needed to 
check several different energies for every direction of H to find which homology class, if any, is associated to 
it. 

Moreover, this means that at every energy different from there is no common point between boundaries 
of different zones, as every zone gets strictly smaller at every change of energy. In the limit for the energy 
that goes to -1 or 1 all zones tend to disappear as above 1 or below -1 the level surface of / is a sphere. 

Let us now examine in more detail the case of energy: it is easy to verify that this surface has 
curvature everywhere negative except in the eight points (±.5, ±.5, ±.5) in which is 0. This means that for 
every direction different from (±1, ±1, ±1) wc will have exactly four critical points, all of saddle type because 
of the topological constraints. 

The analytical expression of the critical points for a generic E G (—1,1) is very complicated but it 
gets much simpler in the most interesting case, namely E = Q. Their expression in cartesian coordinates 
(a, 5) e [0,1]2 is: 



4 Study of the trigonometric function 
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xi{a,b) 
yi{a,b) 

zi(a, 6) 



^ sin ^ (aa (a, b)) 
^ sm~^ {ba{a, b)) 



2-K 



sin ^ {oi{a, b)) , a>b 



{x2,y2,Z2)= - {xi,yi,zi) 

{X3,y3,z3)= {^,^,^) + {xi,yi,zi) 
{xi,yi,Zi)= (1,1,1) - {xx,yx,zi) 

a{a, b) - 



2Va^ + 64 + 1 - a2^2 _ ^2 _ 52 _ (^2 + ^2 + 1) 



a4 + 64 + 1 _ 2a262 _ 2a^ - 2b^ 



It is straightforward to realize that the averaged Euler characteristic xh{c) =< H, 7(c) >, where 7(c) = 
^ WiXi{c) is the sum over all critical points weighted by the Dynnikov index wi equal to the "hamiltonian" 
index of the critical point (as of the 1-forni, see [Dyn97]) times < f,H >, is identically for c = 0. 

This fact is also clear from the relation xh{c) = ^+ — i-^- the averaged Euler characteristic is 

equal to the sum of the height of cylinders of "positive" closed leaves (the ones on which the gradient points 
to the exterior of the loop) minus the height of cylinders of "negative loops" . 

By the symmetry at c = 0, that exchanges "electrons" with "holes" , i.e. positive cylinders with negative 
ones, it is clear that the sum is zero, while it is negative for c > and positive for c < 0. This corresponds to 
the fact that all "ergodic" or "non generic" directions appear just at energy 0, as the nullity of the averaged 
Euler characteristic is a necessary condition for the appearance of these directions. 

That there could be no "ergodic regime" for energies different from was also clear from the fact that 
all energy intervals [ei(i7), e2(-ff)] for which open orbits exist are of the form [— e, e]: "ergodic" directions 
correspond to the case of length zero of this interval, that in this case implies ei = 62 = 0. 

Now let us see what is possible to do "by hand" about stability zones at energy 0. As we have the 
explicit analytical expression for all critical points we can use the following procedure: first of all we make 
sure somehow that a direction (a, b) is "generic" , i.e. it is inside some stability zone, for example looking at 
the plane section generated by the NTC library or by any computer algebra program like Mathematica and 
verifying that just one of the loops is homotopic to (at energy it is enough to examine just one of the 
critical points because of the symmetry). 

Then we choose one critical point, saypi = (a;i,yi,zi), inside the cube [0,1]^ and follow "vertically" the 
cylinder of closed orbits until we reach the second base point. As we observed before, the second base point p 
must be its symmetrical respect to the origin, namely the one we called P4, so in the covering its coordinates 
will be of the form p4 + (I, m, n). Equivalently, going from pi to pi inside the cylinder and coming back to 
Pi through the segment that joins them inside the cube will produce a loop of homology class {l,m,n) in 
T^. As at the boundary of a zone both cylinders have height 0, i.e. the two bases belong to the same leaf, 
it follows that the boundary of any zone is a subset of the curves {< H,pi — P4 — {I, m, n) >= 0}(i_m,n)(£Z3- 

By the topological stability of curves homotopic to 0, this triple of integers depends continuously on 
the magnetic field, so it is locally constant. The number of different triples inside a single stability zone 
determines the number of sides of the zone as shown in figure 4.4. 

The cylinder identified by pi and P4 will disappear either when its height goes to or when it gets 
substituted by a new one: in the first case it means that we reached the boundary of the stability zone. 
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i?= (.53, .268,1) F~ (.5352, .268,1) = (.537, .268, 1) 

Figure 4.2 : the disappearence of a cylinder at the boundary of the stabiUty zone (0,0,1); on the 
left the cylinder has non-zero height, in the central picture the two basis collapse one over the other and in 
the third one the cylinder has disappeared, substituted by a new one relative to the stability zone (1,2,4). 

In the second case it happens that either two different cyUnders colUde and mutually exchange one of 
their bases or a single cylinder collides with itself and the base point is exchanged with one equivalent to it 
but in a different position (see figure 4.3 and tables 1. 1-1.9). 




-0.5 -0.25 0.25 0.5 0.75 1 -0.5 -0.25 0.25 0.5 0.75 1 -0.5 -0.25 0.25 0.5 0.75 1 



iJ = (.22,.23,1) 77= (.23, .23,1) 7? = (.24, .23, 1) 

Figure 4.3 : a change of cylinder inside a stability zone. On the left is shown a critical leaf at the base 
of a cylinder, the critical point is pi — (0.035, 0.463, 0.25). At the opposite base lies the critical point 
p = P4 + (0,0, 1). The middle picture shows what happens at the boundary between the two stability 
zones of cylinders, namely the point pi has a saddle connection with P2- The picture on the right shows the 
base of the new cylinder. At one base still lies the point Pi but at the opposite one now lies ^4 + (1, 1, 0). 

In picture 4.4 we show what happens in case of the zone (2,4,5): there are three different kinds of 
cylinder, labeled by (—3,3,-2), (0,0,1) and (-4,2,-1), so the zone is a triangle divided inside in three 
sub-zones. At the boundary between the first and the second sub-zone the change is determined by the 
appearance of a saddle connection between pi and p2 + (—1, 2, —1), at the boundary between the first and 
the third we have an analogous situation between and^)2 + (— 3, 3, —1) and at the boundary between second 
and third we have instead the appearance of saddle connection of pi with itself, precisely with pi + (2, —1,0). 

In table III and also in the other pictures with smaller resolution it is possible to recognize in many 
stability zones the boundaries between sub-zones in which pi has a saddle connection with itself, as in these 
points the 2-tori filled by open orbits have rank 1 and so these points are not included in the data and the 
stability zone is cut by a segment of straight line. It is easy to check that the same straight line, whose 
equation is la + mb + n = for having a saddle connection with pi + {l,m,n), cuts several (possibly 
infinite) zones. 



4 Study of the trigonometric function 
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Figure 4.4 : the "cylinder" structure of the zone (2, 4, 5). Keeping fixed the critical point of one of 
the bases of the cylinder, say pi, in subzone I in the second base we find the critical point p-i + ( — 3, 3, — 2), 
in subzone II the critical point P4 + (0, 0, — 1) and in subzone III the critical point P4 + ( — 4, 2, — 1). All 
these three cylinders are formed by closed loops that have the same homology class in AIq, namely (2, 4, 5) 
(using coordinates with respect to the natural base in keri*). 

In tables 1. 1-1.9 are shown the three kinds of cylinder corresponding to the three internal subzones and 
the way they transform when the magnetic field direction crosses the internal boundaries. 

These rank-1 2-tori survive longer to energy changes than the rank-2 ones, i.e. we still find them when 
the rest of the zone has disappeared, but they disappear for any generic perturbation of H. A further 
confirniation of the accuracy of our algorithm is given by the perfect agreement between the segment found 
analytically for zone (2, 4, 5) shown above and the one that is possible to see in table III. 

All techniques described above allow us in principle to find analytically all boundaries of stability zones 
and the boundaries of their sub-zones, even though they do nothing to help us finding which homology 
class is associated with them; this quantity of course is anyway easily obtained through our library. The 
main problem is that we did not find any way to put these procedures in any simple algorithm for letting a 
computer do the job, so it has to be done "by hand" . 

Anyway to be able to get this analytical expression does not seem to be crucial in itself: with our NTC 
library we can obtain a good approximation of the interior of any stability zone by sampling the square 
[0, 1]^ with step 1/N in both directions. In that way we will get for every point {m/N, n/N), < m,n < N, 
the homology class of the stability zone it belongs to (if any). It is good though to have such analytical 
expressions as they provide a way to double check the accuracy of our algorithm comparing the interior of 
the zone found with the NTC library with its analytical boundary. 

We initially run our program with resolution N — 100 at energies E ~ 0, — .1, — , 2, — .3, — .5, — .7, — .9 
and found the pictures we show in tables V-X. In table IV are shown the labels associated to the biggest 
zones together with their boundaries found analytically. The boundaries are also drawn in table V to show 
the very good agreement with them of numerical data found using the NTC library. After trying several 
different machines with different operating systems, it turned out that the fastest machines available to us 
were Pentium II Linux machines, so we run all our simulations on them. Every sampling with N = 100 
resolution takes around 12 hours of CPU. 

In table III are shown the data found with the NTC library at a resolution = 1000. The calculation 
explored just the upper triangle 6 > a C [0, 1]^, it run ~ 3 weeks on 5 Linux machines with Pentium II CPUs 
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and found ^ 3 • 10^ distinct stability zones. In tabic III we show just the 1066 zones containing at least 10 
points and then we extend the picture by symmetry to the whole square. In table II we extended the picture 
to the upper half of the sphere by symmetry to show the global pattern of the fractal. 

4.1 Evaluation of the fractal dimension 

One of the most standard procedures to get the fractal dimension of a set is to evaluate its "Box 
Counting" dimension [ASY96,Fal97]. To double check our results we used two different methods to get this 
estimate. 

The first method comes directly from the definition, namely we divide the square in 2^" squares of area 
1/2^" and count how many of them we need to cover the fractal (i.e. the white spots in table III). Below we 
show the data for n = 1, • • • , 10 
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Figure 4.1 : plot of the log in base 2 of number of squares needed to cover the fractal with squares 
of size 2^" versus the size scale n. 



After we discard the last two terms, that probably wc cannot evaluate well enough because of the 
finiteness of our resolution, we find that the slope that minimizes the rms in a linear fitting of the above 
plot is d~ 1.78. 

The other method wc used is the following: after having fixed a number r > 1 we count for every n 
how many zones have area between r~" and r~"~^. Let us call this number A^^. Then as n ^ oo the ratio 
between log^{Nn) and n converges, for fractals for which that dimension is well defined, to the Box Counting 
dimension divided by the dimension of the ambient space [Fal97]. 

The picture below shows the plot in case r = 2: 
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Figure 4.2 : plot of the log in base 2 of number of number of zones of area between 2 " and 2 
versus the size scale 71. 



-n-1 



In this case the global behavior is much less linear but is clear that the first points have no real meaning 
because there the scale is still too big and it is safe also to discard the last ones as there we are probably at 
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a scale too small for the resolution of our picture. After discarding these boundary points we are left with a 
plot which can be well approximated by a linear function with slope a ~ .91. This suggests that d ~ 1.82, 
which is in very good agreement with the previous estimate for the fractal Box Counting dimension of the 
set of "ergodic" directions. 

5 Study of the piecewise quadratic function 

Using degree-2 polynomials we can build a function that has the same properties of the previous one but 
is much easier to deal with analytically. The function will not be globally smooth as the second derivatives 
will not glue smoothly, but still it will be globally and piecewise quadratic. 



In particular it is possible to build a piecewise polynomial function of degree 2 that allows to evaluate the 
espression of all critical points at every energy, so that we will be able in principle to verify the agreement of 
our algorithm with every zone at energies different from zero and to find analytical expressions for topological 
quantities that depends on them like the averaged Euler characteristic. 

The function we used is the following: 



where [x] is the fractional part of x for a; > and F is extended to (— oo,0) by F{—x) = —F{x). 

Its level sets are very similar to the ones of previous function. Below we show a picture of the level 
No = /~^(0), that has the same peculiarity of the level set Mq studied in previous section. 

As before, this function in the range of energies (—1, 1) gives rise to genus-3 surfaces embedded in T'^ 
with rank 3, so just 4 saddles of the foliation contribute to the topology of our system. All other saddles (if 
any) will be linked to a center and hence will be homotopic to in the surface and easily eliminated from 
the surface through a homotopy naturally generated by the center itself. 

The analytical expression for the critical points for E' e [—1, 0] are the following: 




Figure 5.1 : the surface A^o = / ^(0). The three basic cycles non homotopic to in T' 



are shown. 
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xi{a,b,E)= < 



yi{a,b,E)= < 



Zi{a,b,E)= < 



X2{a, b, E)- 



y2{a,b,E)= < 



Z2{a,b,E)-- 
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, b''-a^>E 
, b"^ -a^ <E 

b^-a^>E 

b^-a^<E 

b^-a^>E 

a^-b^<E 

a^-b^<E 
a^-b^>E 

, -b"^ <E 
a^-b^>E 

a^-b^ <E 
a^-b^>E 



X4= 1 — Xi 



2/4= 1 - 2/1 



Zi = 1 — Zi 



X3= 1- X2 



ys = 1 - 2/2 



-23 = 1 - -22 



The expression of boundaries of all zones in this case is very simple: for example the boundary of the 
zone labeled by (0, 0, 1) at energy E is the union of the segments of ellipse (1 + E)a'^ + (3 — -E')62 = 1 + for 
b> a and (1 + E)b'^ + (3 — E)a'^ = 1 + E for b < a, and the boundary of the zone corresponding to (1, 1,1) is 
the union of the segments 8a - (3 - £)a2 - (1 + i;)62 = 3 - ^ for 6 > a and 86 - (3 - £;)62 - (1 + E)a2 = 3 _ 
for b < a. 

Using the same triples of integers used for boundaries in table IV we have been able to find with a 

few modification the corresponding zones for this function. As shown in table XI to this zones corresponds 
exactly the same homology classes of the previous picture, as we expected given the similarity between the 
two functions. 

We analized numerically the stability zones in the square [0, 1]2 for the same energies, finding the data 
reported in tables XII- XVIII. At every energy we included in the picture also the boundary of a few zones 
to show the very good agreement of numerical data with the analytical results. 

Using the data found at resolution 1000 (table XI) we evaluated again fractal dimension of the set of 
crgodic directions with the two methods used for the trigonometric case, finding very similar results: the 
Box Counting method gives us an estimate of d ~ 1.77 and from the growth rate of sizes of stability zones we 
get d/2 ~ .9. Therefore the two different estimates are in very good agreement also in this case and suggest 
a fractal dimension around d = 1.8. 
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Figure 5.2 : plot of the log in base 2 of number of squares needed to cover the fractal with squares 
of size 2^" versus the size scale Tl. 
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Figure 5.3 : plot of the log in base 2 of number of number of zones of area between 2"" and 2 

versus the size scale n. 



6 Conclusions 

We produced a C++ library that implements all functions needed to analyze numerically the topological 

behaviour of foliations induced on a periodic surface of genus 3 by a constant 1-form. This problem is 
equivalent in physics to the behaviour of magnetoresistance under a strong magnetic field. 

We checked numerically our code on two "toy functions" that produce genus 3 surfaces embedded in 
with rank 3 and verified its correctness comparing numerical data with the analytical data that was possible 
to get for the two simple functions chosen, finding a very good agreement between the two. 

Our next move will be to apply this machinery to concrete Fermi surfaces of normal metals, that have 
genus 4 in the easiest cases. 
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1. = (.37, .73, 1) 



2. Hc^ (.37, .742, 1) 



3. H= (.37, .743, 1) 




4. H= (.41, .77, 1) 5. H = (.385, .77, 1) 6. H = (.384, .77, 1) 




7. = (.42, .76, 1) 8. if ~ (.42, .761, 1) 9. = (.42, .77, 1) 

Table I : to illustrate the phenomenon of the ehange of cylinder type inside a stability zone we show what 
happens in case of the zone (2, 4, 5) shown in figure 4.4 . 

1-3: in the first row we move the direction of H from subzone Ito subzone III. In the central picture we reach 

the boundary between the subzones, in which point pi has a saddle connection with point P2 + ( — 3, 3, — 2). 
4-6; in the second row we move H from subzone Illto subzone II. This time on boundary pi has a saddle 
connection with a copy of itself separated by a 1-rational vector (2, —1, 0). 

7-9: in last row we move H from subzone Ito subzone II. At the boundary pi has a saddle connection with 

P2 + (1,-2,1). 



Table II : the fractal picture in the square [0, 1]^ obtained at a resolution N = 10"^. Of the ~ 3 • 10^ 
zones found just the ones with at least 10 points ('^ 1000) are shown. The square has been obtained just 
symmetrizing the triangular picture obtained. To get this picture we used 5 Linux machines with Pentium II 
CPUs for 3 weeks. It is possible to get the homology class corresponding to the biggest zones comparing 
this picture with next one and with the table included in next page. From these data has been extrapolated 
a fractal dimension of C? ~ 1.77 for the set of "ergodic" directions. 



Table III : the fractal picture in the square [0, 1]^ obtained at a resolution N = 10^. Of the ~ 3 • 
zones found just the ones with at least 10 points ('^ 1000) are shown. The square has been obtained just 
symmetrizing the triangular picture obtained. To get this picture we used 5 Linux machines with Pentium II 
CPUs for 3 weeks. It is possible to get the homology class corresponding to the biggest zones comparing 
this picture with next one and with the table included in next page. From these data has been extrapolated 
a fractal dimension of d — 1.77 for the set of "ergodic" directions. 
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Table IV : boundaries of a few stability zones together with their homology class. All these boundaries 
have been obtained with the technique explained in section 4, while the associated homology class has been 
obtained numerically, except in the trivial case of (0,0,1). Below are listed the zones with biggest size and 
their area, form the data found at iV = 1000. 



Horn Class 


Area 


Horn Class 


Area 


Horn Class 


Area 


(0,0,1) 


(2.83 ± .02)10-^ 


(2,4,5) 


(8.6 ± .6)10-^ 


(1,6,6) 


(2.0 ± .1)10-3 


(1,1,1) 


(2.03 ± .01)10-^ 


(1,4,4) 


(8.3 ± .3)10-^ 


(4,5,8) 


(2.0 ± .4)10-3 


(1,2,2) 


(8.2 ± .2)10"^ 


(1,2,4) 


(6.2 ± .5)10-^ 


(5,8,10) 


(1.9 ± .4)10-3 


(0,1,2) 


(5.1 ± .1)10-^ 


(3,4,6) 


(4.7 ± .5)10-^ 


(4,6,9) 


(1.8 ± .3)10-3 


(1,3.3) 


(2.1 ± .1)10-2 


(1,5,5) 


(4.1 ± .2)10-3 


(1,6,10) 


(1.7± .1)10-3 


(2,3,4) 


(1.7± .1)10-2 


(2,5,8) 


(4.1 ± .4)10-3 


(5,9,11) 


(1.6 ± .2)10-3 


(1,3,5) 


(9.6 ± .5)10-^ 


(4,7,8) 


(3.0 ± .3)10-3 


(4,6,7) 


(1.5 ± .2)10-3 


(1,4,6) 


(9.6 ± .5)10-^ 


(0,3,4) 


(2.9 ± .4)10-3 


(2,3,6) 


(1.5 ± .4)10-3 


(0,2,3) 


(9.0 ± .6)10-^ 


(3,5,7) 


(2.7 ± .3)10-3 


(3,5,9) 


(1.5 ± .4)10-3 
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Table V : map of the stability zones at energy in the square [0, l]'^ sampled at a resolution N — 100. 
Nearly 700 zones are found at this resolution; in the above picture we show just the 74 that contain at least 
5 points. The boundary found analytically is also shown for a few zones to show the perfect agreement with 
the numerical results. 
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Table VI 

- 100. 

point. 



: map of the stability zones at energy E = —.1 in the square [0, 1]^ sampled at a resolution 
Just 48 zones are left at this energy, and here we plotted just the 34 ones with more than 1 
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Table VII : map of the stability zones at energy E — — .2 in the square [0, l]'^ sampled at a resolution 
N — 100. Just 12 zones are left, and we show all of them in the picture above. 
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Table VIII : map of the stability zones at energy E = — .3 in the square [0, 1]'^ sampled at a resolution 
N ~ 100. All 8 zones found are shown. 
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Table IX : map of the stability zones at energy E — — .5 in the square [0, 1]^ sampled at a resolution 
N ~ 100. At this energy just the four zones shown arc left. 
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Table X : map of the stability zones at energy E = — .9 in the square [0, 1]^ sampled at a resolution 
N — 100. Just the two biggest zones are now visible. 



Table XI : the fractal picture for the piecewise quadratic function in the square [0, 1]^ obtained at a 
resolution TV = 10^. Of the ~ 3 • 10^ zones found just the ones with at least 10 points {'^ 1000) are 
shown. The square has been obtained just symmetrizing the triangular picture obtained. 
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Table XIII : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at _E/ — with resolution N — 100. Of the 708 zones found, just the 74 with more than 5 points are 
shown. The boundary found analytically as explained in section 5.1 is also shown for a few zones to show 
the perfect agreement with the numerical results. They are very close to the boundaries of trigonometric 
function shown in table IV and the homology zones that labels them are exactly the same than in the 
trgonometric case. 




Table XIV : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at E = —.1 with resolution N — 100. We show here all 42 zones found together with the boundaries of 
the biggest ones. 




Table XV : map of the stability zones for the piecewise quadratic function in the square [0, 1]"^ sampled 
at E =^ —.2 with resolution N — 100. All 14 zones found are shown together with boundaries of the 
biggest ones. 
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Table XVI : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at E ^ —.3 with resolution N — 100. All 10 zones found are shown together with boundaries of the 
biggest ones. 
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Table XVII : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at — — .5 with resolution N = 100. All 4 zones found are shown together with their boundaries. 
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Table XVIII : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at — — .7 with resolution N = 100. Just the two main zones survive at this energy. 




Table XIX : map of the stability zones for the piecewise quadratic function in the square [0, l]'^ sampled 
at = — .9 with resolution N = 100. Just the two main zones survive at this energy. 



